Numerical Study of Flow Past a Circular Cylinder 
Using RANS, Hybrid RANS/LES 
and PANS Formulations 


Alaa Elmiligui* 

Analytical Services & Materials, Inc. 
107 Research Drive, Hampton, VA 23666 

Khaled S. Abdol-Hamid| 

NASA Langley Research Center 
Hampton, VA 23681 

Steven J. Massey J 
Eagle Aeronautics, Inc. 

13 W. Mercury Blvd., Hampton, 23669 

S. Paul Pao § 

NASA Langley Research Center 
Hampton, VA 23681 


Abstract 

Two multiscale type turbulence models are implemented in the PAB3D solver. The 
models are based on modifying the Reynolds Averaged Navier-Stokes (RANS) equations. 
The first scheme is a hybrid RANS/LES model utilizing the two-equation (kg) model with a 
RANS/LES transition function dependent on grid spacing and the computed turbulence 
length scale. The second scheme is a modified version of the partially averaged Navier- 
Stokes (PANS) model, where the unresolved kinetic energy parameter (f k ) is allowed to vary 
as a function of grid spacing and the turbulence length scale. Solutions from these models 
are compared to RANS results and experimental data for a stationary and rotating cylinder. 
The parameter f k varies between zero and one and has the characteristic to be equal to one 
in the viscous sub layer, and when the RANS turbulent viscosity becomes smaller than the 
LES viscosity. The formulation, usage methodology, and validation example are presented to 
demonstrate the enhancement of PAB3D T s time-accurate and turbulence modeling 
capabilities. The models are compared to RANS results and experimental data for turbulent 
separated flows (TS) and laminar separated flows (LS) around stationary and rotating 
cylinders. For a stationary cylinder, the TS case is accurately simulated using the general 
two-equation kg turbulence model (eddy-viscosity model). PAB3D accurately predicts the 
drag coefficient (C D ), lift coefficient (C L ) and the Strouhal number (St). The LS case was a 
challenge for the RANS computation with an eddy-viscosity turbulence model. The 
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RANS/LES and PANS performed well and showed marked improvements over the RANS 
solution. The modified PANS model was the most accurate. For the rotating cylinder, the 
spin ratio varied from zero to one, and the PANS results were in good agreement with 
published experimental data. RANS/LES and PANS capture both temporal and spatial 
fluctuations and produce large-scale structures that do not occur in the RANS simulation. 
The current results show promise for the capability of PANS in simulating unsteady and 
complex flow phenomena. 


I. Introduction 

The limited capability of the Reynolds Averaged Navier-Stokes (RANS) approach combined with eddy- viscosity 
turbulence models to simulate unsteady and complex flows has been well known for some time. The RANS assumes 
that most of the energy is modeled through the turbulence transport equations and is resolved in the grid. RANS also 
over predicts the eddy viscosity, which results in excessive damping of unsteady motion. Consequently, the eddy 
viscosity of the unresolved scales attains unphysicaly large values suppressing most temporal and spatial 
fluctuations in the resolved flow field. One of the approaches to overcome this problem is to provide the required 
mechanism to resolve the largest scales of motion. Among several methods, the Detached Eddy Simulations (DES) 
[1], hybrid Large Eddy Simulation (LES) [2-3], Limited Numerical Scheme (LNS) [4] and Partial Averaged Navier- 
Stokes (PANS) [5] are capable of providing the needed mechanism to satisfy this requirement. One of the major 
deficiencies associated with the heretofore published use of hybrid schemes is that there is no clear identification of 
the different flow regions. These regions need to be clearly defined as RANS regions and hybrid regions in order to 
achieve complete simulation, independent of grid resolution. Several researchers observed that in most cases, using 
hybrid methods, the use of fine grid might result in incorrect simulations. Abdol-Hamid and Girimaji [6] explored 
new approach to improve the accuracy and robustness of creating a simulation of an unsteady flow field based on 
the work by [6]. They accomplished that through the development and implementation of a novel two-stage 
procedure to efficiently estimate the level of scale resolution possible for a given flow on a given grid for PANS and 
other hybrid models. 


PAB3D is a structured, multiblock, parallel, implicit, finite- volume solver of the three-dimensional RANS 
equations, and advanced turbulence models are available in the code. PAB3D is widely used for internal and 
external flow applications by NASA and by the US aerospace industry. Investigations in the area of unsteady flow 
control for propulsion applications have led to an increased interest in upgrading PAB3D’s [7-8] time-accurate 
capabilities. The current version of the PAB3D code has a second-order time accuracy algorithms scheme by 
employing either physical time with sub-iteration or dual time sub-iteration [8]. 

In an attempt to increase the flow physics fidelity and the accuracy of PAB3D code, hybrid turbulence model 
RANS/LES [2-3] has been added. An alternate new feature to PAB3D is the addition of the Partially Averaged 
Navier-Stokes (PANS) method, which was suggested by Girimaji et al. [5]. The PANS model was developed to 
overcome the grid dependency associated with the customary implementation of the hybrid RANS/LES method. The 
addition of improved algorithms for second-order time accuracy, sub-iteration schemes, hybrid turbulence modeling, 
and moving boundary conditions provide key modernizing enhancements to the code. The primary objective of this 
paper is to assess the performance of these newly implemented approaches in a production CFD code. 

The organization of the paper is as follows: The governing equations of the RANS, RANS/LES and PANS 
formulation will be presented and discussed in detail. Computational results from RANS, RANS/LES and PANS for 
a flow past a stationary and a rotating cylinder will be presented, and compared to experimental data. Flow around a 
cylinder is considered as the test case for the hybrid turbulence model [9-13], because it is a basic engineering 
problem and is inherently unsteady. 


II. Approach 

The governing equations of the RANS formulation include the conservation equations for mass, momentum, 
energy, and the equation of state. In the present study, the perfect gas law is chosen to represent the air properties, 
and the eddy viscosity concept is used to model the Reynolds stresses. The mass, momentum, and energy 
conservation equations of the RANS equations can be written in a conservative form as follows: 
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The boundary conditions for s and k at the wall are: 
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and k wa n = 0. The turbulent stress components are formulated as: 
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For the purpose of this paper, we will define RANS turbulent viscosity as 

k 2 


o MNS = f M pC, 


(4) 


A. Hybrid RANS/LES 

Nichols and Nelson [2] give an example of a hybrid RANS/LES turbulence model. This method was 
implemented in conjunction with Menter’s SST two-equation turbulence model and is termed a multiscale (MS) 
model. In the present paper, the hybrid model is used with the two-equation model described in equations 2 and 3. 
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The turbulent length scale, used in this implementation, is defined as 


l t =max(6.*Jv t /Q ,k 312 Is) (5) 

The sub grid turbulent kinetic energy is defined as 

k LES =f d k (6) 

The damping function is defined as 

f d = {1 . + tanh[2;r(A - 0.5)]}/2 (7) 


where, 
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X is the unresolved characteristic ratio, and 

A = max(A x ,A y ,A z ) (11) 


The eddy viscosity is then calculated from: 


o t =f d or s Hi.-f d )or 

(10) 

v LES =min(uf ANS ,0M4A^) 

(12) 


Note that this hybrid model allows the transition from RANS to LES as a function of the local grid spacing and the 
local turbulent length scale predicted by the RANS model rather than as a function of the grid spacing alone. This 
allows the model to detect whether it can resolve the turbulent scales present on the existing grid before its transition 
over to the LES mode. 

B. PANS Approach 

In its original form, PANS [5] replaces the two-equation turbulence model by solving for the unresolved kinetic 
energy k u and the dissipation S u . The k u equation is identical to the original k equation. The dissipation equation has 
only one major change through: 

C £ 2 - fk(C £2 — Qi)+ C £l (13) 

In the present work, we introduce an attempt to use a variable f k instead of a constant value. We use equation (7) 
to compute f k as: 

f k = { 1 . + tanh[2;r(A - 0.5)]}/2. (14) 

In this case, turbulent length scale is defined as: 

K = k l' 2 A = 1 + ^ 4 / 3 ’ and A = T 

The function in equation (14) has the characteristic to be equal to 1.0 in the viscous sub layer, as the unresolved 
characteristic ratio tends to be of very small value. Also, the value of this function is restricted to 1, in case the 
RANS turbulent viscosity becomes smaller than the LES viscosity (12). 
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C. Boundary Conditions 

In the present study, a characteristic Riemann invariant type boundary condition was used to model the far field 
boundaries, while a periodic boundary condition was imposed in the longitudinal direction of the cylinder. A no-slip 
velocity boundary condition was used at the cylinder surface so that the fluid would neither slip nor penetrate. For 
the case of the stationary cylinder, there is no relative motion between the cylinder surface and the fluid, and the 
velocity components were set to zero. For the rotating cylinder, the velocity at the cylinder surface was set to be 

equal to COX r where CO is the cylinder angular velocity and r is the position vector connecting the cylinder axis of 
rotation to the cylinder surface as shown in Fig. la. The interface to the new rotating/spinning boundary condition 
resides in the user.cont file. The user has to specify the angular velocity (rad/sec), the axis of rotation, and a point on 

the axis. Vector r is the vector with the minimum distance between the rotating surfaces to the axis of rotation, and 
is defined as: 

d = |(*2 -*l) X (*i -X 0 )| 

\x 2 - x x | 

where x\ &x 2 are points on the axis of rotation while x 0 is a point on the surface of cylinder. 


III. Results and Discussions 

Flow past a stationary cylinder as well as a rotating circular cylinder, was computed to verify the time accuracy 
of the code and of the relative advantage of the hybrid turbulence models RANS/LES and PANS. The two- 
dimensional grid consisted of 32,256 cells and 24 blocks, and extended 15 diameters into the far field (Fig. lb). The 
three-dimensional grid was the same as the two-dimensional grid with the addition of 40 planes, which covered two- 
cylinder diameter. The same grid was used for all runs, which gave a first grid y+ range of approximately .2 to 2.0. 
The diameter of the cylinder D was 40 mm wide at Re = 50,000 and the Mach number for all cylinder cases was set 
at a value of M = 0.3. A non-dimensional time step of 0.015 (based on free stream speed and the diameter of the 
cylinder) was used for all cases. Based on the Strouhal number range and the time step used for the present cases, 
approximately 350 data points in time per cycle of shedding was sampled. Four sub-iterations were used to reduce 
the error. Each of the 2D simulations required approximately 3 hours using 24 (2.8 GHz P-4) computers. The 3D 
cases each required approximately 48 hours using the same set of 24 computers. 

A. Stationary Circular Cylinder in Cross-flow: 

The shedding frequencies are determined from the lift coefficient (C L ) fluctuation as it varies with time (Fig. 2). 
The C L is obtained via the internal force and moment integration algorithm, Post [14]. The strategy for this case was 
to first run the simulation at a very coarse grid level ( 1 /4 th in each direction) for 10,000 iterations to trigger its 
asymmetric vortex shedding instability, and then refine the grid and run the solution for an additional 20,000 
iterations. The solutions were averaged over the last 15,000 iterations (approximately 50 shedding cycles). The onset 
of asymmetric vortex shedding is seen to occur just after the first 60 time unit, and the switch to fine grid is seen to 
coincide with an increase in amplitude. It was observed that approximately 4 sub-iterations per physical time step 
produced the optimal convergence per iteration. However, the physics of the specific problem will dictate the sub 
iteration number for other cases. In the present results, four sub iterations typically reduced the residual by three 
orders of magnitude at that time level, with no improvement using more iteration. The results were compared with 
the results using up to 20 sub-iterations, with no substantial difference in the final results. 

The Strouhal number (St) was captured with the use of the Power Spectral Density (PSD) of C L as shown in Fig. 
3. To verify the capability of PAB3D for simulating unsteady flow problems, the Turbulent Separated (TS) case was 
used, which simulates the experimental flow condition in which the boundary layer is tripped well ahead of 
separation. Similar to what was found by other researchers [15-17], we achieved this objective numerically by 
choosing a free stream turbulence level high enough to cause natural transition (5 times laminar viscosity). We 
performed 2D-RANS computations with PAB3D using two-equation ks at the Reynolds number of 140,000. This 
flow condition matched the conditions used by Travin et al. [15], Vatsa and Singer [16] and Hansen and Forsythe 
[17]. The results from the work by Travin et al. [15] show that there are only small differences observed between the 
2D RANS, 3D RANS and 3D-DES results. For this reason, we will not present any 3D results for this case. The 
time-averaged surface pressure coefficients (C P ) resulting from the 2D computations were compared with the 
experimental data of Roshko [10] in Fig. 4. Table 1 compares the present PAB3D results using the 2D kc turbulence 
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model with the 3D-DES studies done by Travin et al. [15] and Hansen and Forsythe [16] and the experiment of 
Roshko [10]. Based on the results listed in Table 1 and the comparisons of Fig. 4, it was concluded that the current 
RANS compared well with Roshko’ s experimental data and with other CFD simulations. 


Table 1. Re=140,000 TS Time-Averaged Results 


Method 

Cd 

Cpb 

St 

3D-DES Travin et al. [15] 

0.57 

0.65 

0.3 

Hansen and Forsythe [17] 

0.59 

0.72 

0.29 

PAB3D ks 

0.62 

0.68 

0.27 

Experiment 

0.62-0.74 

0.5-0. 9 

0.27 


The laminar separated (FS) case was chosen to evaluate the implementations of both Hybrid RANS/FES and 
PANS formulations into the PAB3D code, and to investigate their capability to simulate such flow. This case is 
compared with the experimental data of Nordberg [18] (Re=3000 and 8000) and Cantwell and Coles [19] 
(Re=140,000). Fig. 5 shows the PSD variation with St from RANS and PANS results. It is clear that RANS 
produced a peak at only one major Strouhal number and no significant spectral energy at higher Strouhal numbers. 
This result indicates that the 3D RANS simulation has a two-dimensional flow character. The vorticity plot shown in 
Fig. 6 supports this conclusion as it shows the flow with no change in the Y-direction. On the other hand, the PANS 
resulted in one major Strouhal number and several minor ones, indicating that more scales were resolved in this 
simulation. Similar observations were made for the RANS/FES simulation. Fig. 7 shows the three-dimensional 
character of the flow produced as a result of the PANS formulation. Table 2 shows the comparison between the 
recent PAB3D simulations, Travin et al. [15] DES solutions, and the experimental data. The PANS formulation in 
PAB3D produces the closest results as compared with the experimental data. The 2D PANS and FES results largely 
differ from the experimental data, but are similar in character to the 2D-DES simulations of Travin et al. [17]. 

Table 2. Re=50,000 LS Time- Averaged Results 


Method 

C D 

-Cpb 

St 

3D-DES Travin et al. [15] 

1.27 

1.28 

0.21 

2D-DES Travin et al. [15] 

1.77 

2.05 

0.14/0.20 

2D PAB3D RANS/LES 

1.69 

2.05 

0.21 

2D PAB3D PANS 

1.67 

2.1 

0.22 

3D PAB3D RANS/FES 

1.0 

0.9 

0.22 

3D PAB3D PANS 

1.1 

1.03 

0.21 

2D/3D PAB3D RANS 

1.08 

1.03 

0.23 

Experiment 

0.98-1.25 

0. 9-1.2 

0.18-0.21 


Figure 8 shows the variations of C L with respect to time for the PANS formulation. The 3D PANS solution 
displayed large and random variations, as a result of strong modulation of the vortex shedding. These observations 
are similar to the ones reported by Travin et al. [15] and Vatsa and Singer [16], using 3D-DES formulation. The 
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time-averaged surface pressure distributions from the 3-D simulations are compared with experimental data of 
Nordberg [18] and Cantwell and Coles [19] in Fig. 9. In general, the hybrid RANS/LES and PANS formulations 
produced comparable results to the 3D-DES solutions presented in Refs. [15 & 16]. 

Figure 10 shows the unresolved turbulent kinetic energy resulting from the RANS and PANS formulations. The 
levels from the RANS simulation were ten times those produced by the PANS formulation. More energy was 
resolved using the PANS formulation as compared with that produced by the RANS simulations. Fig. 11 shows the 
time-averaged f k from the PANS simulations. The f k was as low as 0.3 in the region of interest where more energy 
was resolved through the use of the PANS formulation. The smaller the f k value is, the more LES-like simulations 
are produced. The closer f k is to a value of 1, the more RANS-like the simulation results. This occurs in the viscous 
sub-layer and far field regions, as expected, where the grid resolution is enough for the use of RANS. 

Rotating Circular Cylinder in Cross-flow: 

The computational grid described earlier and shown Fig. lb was used to compute the flow field around a rotating 
cylinder. The flow field has been computed both in 2D and 3D modes using the standard kc turbulence model and 
the hybrid turbulence models RANS/LES and PANS. Three main parameters govern the flow around rotating 
cylinders: the Reynolds number, the ratio of peripheral to free stream velocity, and the cylinder aspect ratio. 

In the present calculation, the rotational speed of the cylinder varied from 0 to 10000 rpm; the free stream Mach 
number was 0.3 and the Re=50,000. The same solution strategy used to compute the flow field around the stationary 
cylinder was used for the rotating cylinder cases. The variation in the time-averaged C L with respect to the spin ratio 
(SR) is shown in Fig. 12a. SR is defined as the ratio between the peripheral velocity (u) of the cylinder to the free 
stream velocity (U). The rotation of the cylinder generated lift, and C L increased with increasing SR. Comparison 
between the present calculations and the results of Aoki et al. [12] and Tokumaru et al. [13] showed good 
agreement. Simulations using 3D Hybrid RANS/PANS yielded results that are closer to the experimental data [12]. 
Variation of the C D with respect to the SR is shown in Fig. 12b. For the range of speeds computed in the present 
study, the coefficient of drag tends to decrease with increasing SR. The St number for the rotating cylinder was 
captured using the PSD of C L as explained earlier. The range of St for the rotating cylinder varied between the 
values of 0.23 and 0.28. The variation of St with SR is shown in Fig. 12c. Increasing the SR tends to increase the St 
number; however, there was no significant change in the St between the 2D and 3D calculations. 

Figure 13 shows the time-averaged surface pressure distribution as a result of the 3D PANS simulation for the 
SR of 0.0, 0.3, 0.6, and 0.9 respectively. For the stationary cylinder (SR = 0.0), Cp distribution is symmetrical 
around the cylinder. As the SR increased, the pressure distribution became asymmetric due to the spinning, and the 
fact that the net turning of the flow produced lift. Corresponding pressure contours are shown in Fig. 14. The free 
stream flow over the top of the cylinder follows the induced flow from the spinning, while the free stream below the 
cylinder is opposed by the induced flow. This result in an accelerated flow on the top half of the cylinder. 

IV.Concluding Remarks 

A hybrid turbulence model Reynolds Averaged Navier- Stokes/Large Eddy Simulation (RANS/LES) and Partial 
Averaged Navier-Stokes (PANS) have been added to PAB3D code. The new capabilities improve the accuracy of 
simulating an unsteady flow field. In this paper, a new approach to prescribe the unresolved kinetic energy 
parameter (f k ) function was proposed. The parameter f k is a function of length scale and grid size, which represents a 
characteristic length scale. Some of the drawbacks of such a function are that it varies with time and space and that it 
could be affected by grid resolution. The PANS approach is much simpler to implement in most CFD codes than 
other approaches such as DES and the RANS/LES. The preliminary results provide a preview of the potential 
capability of PANS in simulating unsteady and complex flow phenomena. The PAB3D code utilizing 3D PANS 
captured both temporal and spatial fluctuations, and solutions compared well with the experimental data. 

The turbulent separated (TS) and laminar separated (LS) flows were simulated for a stationary cylinder. The TS 
case was accurately simulated using a general two-equation kc turbulence model (eddy- viscosity model) without any 
enhancement. PAB3D accurately predicted the drag coefficient (C D ) and the Strouhal number (St). The LS case 
posted a great challenge for the eddy- viscosity turbulence models. The RANS/LES and the PANS approaches in the 
CFD code provided much better predictions when compared to measured data. We have observed much better 
predictions of C D , Strouhal number, and the surface pressure distribution than the results with the kc turbulence 
model. 

In the case of rotating cylinder, the present calculations indicate that with an increase of the rotational speed, C L 
increases while C D decreases. The PANS provided the best comparison with the experimental data. The computed St 
number is in the range of 0.22 to 0.28, which is in good agreement with previously published results. 
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Future work will involve explorations of new approach to improve the accuracy and robustness of creating a 
simulation of an unsteady flow field. This can be accomplishes through the implementation of a novel two-stage 
procedure to efficiently estimate the level of scale resolution possible for a given flow on a given grid for PANS and 
other hybrid models. 
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Figure 2. The Lift Coefficient (C L ) Fluctuations with Time using 2D PANS Formulation. 
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Coefficient of Pressure 



Figure 3. Power Spectral Density vs Strouhal Number using ke RANS 
Formulation for Re=140,000. 



Figure 4. Comparison between Coefficient of Pressure on the Cylinder Surface 
using 2D ks RANS and Experimental Data of Roshko’s (TS Case). 
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and ks RANS Formulations (LS Case). 
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Figure 6. Three-Dimensional Vorticity Magnitude Results 
from RANS Simulation (LS Case) 
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Figure 7. Three-Dimensional Vorticity Magnitude Results from PANS 
Formulation (LS Case). 



Figure 8. Lift Coefficient (C L ) Fluctuations with Time Results 
from 3D PANS Formulation (LS Case). 
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Figure 9a. Coefficient of Pressure on the Cylinder Surface Results from 3D ks RANS, 
RANS/LES and PANS Formulation at Y=1.0 (LS Case). 
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Figure 9b. Coefficient of Pressure on the Cylinder Surface Results 
from 3D PANS Formulation Compared with Experimental Data (LS Case). 
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Figure 10a. Unresolved Turbulent Kinetic Energy, k u (m 2 /sec 2 ), 
Contours at Y=1.0 Produced from RANS Simulation (LS Case). 



Figure 10b. Unresolved Turbulent Kinetic Energy, k u (m 2 /sec 2 ), Contours at Y=1.0 
Produced from PANS Simulation (LS Case). 
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Figure 11. Averaged f k Contours at Y=1.0 from PANS Simulation (LS Case). 
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Figure 12a. Lift Coefficient for Various Spin Ratios 



Figure 12b. Drag Coefficient for Various Spin Ratios 



Figure 12c. Computed Strouhal Number 
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Figure 14. Averaged Pressure contours for 3D PANS for SR=0.0, 0.3, 0.6 and 0.9 
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